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ABSTRACT 

Two  numerical  methods  of  determining  optimum 
paths  in  the  problem  of  Bolza  are  presented.  The  basic 
theory  of  the  methods  is  outlined  and  their  essential 
characteristics  are  demonstrated  by  several  specific 
applications  to  ship  routing. 

The  first  method  is  characterized  as  one  of  differ- 
ential corrections  employing  methods  developed  by  Professor 
Faulkner  of  the  U,  S.  Naval  Postgraduate  School.   The 
second  method  is  termed  a  method  of  gradients  employing 
calculation  routines  developed  by  Professor  Bryson  of 
Harvard.   On  the  basis  of  the  applications  comparisons  are 
made  in  the  areas  of  ease  of  applicatior.  and  speed  in 
producing  results. 

I  wish  to  express  my  appreciation  for  the  guidance 
and  encouragement  of  Professor  Faulkner  and  for  the 
programming  of  Mary  Haynes  of  the  Postgraduate  School 
Computer  Center. 
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INTRODUCTION 

The  problem  of  obtaining  maximum  or  minimum  performance 
in  a  variety  of  controllable  situations  is  undoubtedly  as  old 
as  time.   Usually  the  degree  of  success  is  determined  by  the 
efficiency  with  which  some  vital,  limited  quantity  is  ex- 
pended . 

Often  the  best  solution  cannot  be  defined,  except  as  it 
develops  by  experience,  even  though  the  problem  can  be  de- 
fined and  stated  mathematically.  The  magnitude  of  the 
analysis  involved  limits  the  feasibility  of  determining  the 
solution.  Modern  high  speed  computing  machinery  has,  for 
many  problems,  removed  the  latter  limitation.   As  a  result 
renewed  interest  has  been  shown  and  much  progress  made  in 
the  development  and  application  of  control  theory  to  such 
problems  of  optimum  control. 

This  paper  presents  an  outline  of  two  practical  nu- 
merical methods  of  solution  of  some  general  problems  of  this 
type.   The  class  of  problems  is  more  specifically  termed  the 
problem  of  Bolza  [l].   For  purposes  of  identification  in  this 
paper  the  two  methods  will  be  called  the  'differential 
method'  and  the 'method  of  the  fundamental  lemma'  which  will 
henceforth  be  abbreviated  the  D-method  and  the  FL-method 
respectively. 

The  D-method  is  based  on  the  formulas  for  differentials 
due  to  Bliss  [2].  The  numerical  scheme  for  determining  the  con- 
stants of  the  solution  by  a  Newton  iteration  is  the  work  of 


Faulkner  [3], [4],   His  research  papers  and  related  studies 
are  the  source  of  both  the  theory  and  the  numerical  data 
presented  for  the  D-method. 

In  the  literature  the  FL-method  is  also  referred  to  as 
the  'gradient1  and  the  'steepest  ascent'  method.   The  calcu- 
lational  procedure  is  similar  to  the  proof  of  the  fundamen- 
tal lemma  and  is  apparently  conceptually  due  to  Courant  [5]. 
The  numerical  scheme,  which  makes  possible  its  application 
to  the  class  of  problems  being  considered,  is  the  work  of 
Bryson  and  Denham  [6],   A  similar  method  is  presented  by 
Kelly  [7],   The  theory  and  basic  procedure  for  application  of 
the  FL-method,  as  reported  in  this  paper,  is  that  of  Bryson 
and  Denham.   The  contribution  of  the  writer  lies  in  the  de- 
velopment of  specific  convergence  schemes,  left  undefined  by 
the  authors  [6],  which  make  possible  the  programming  and  so- 
lution of  two  relatively  simple  ship  routing  problems  thereby 
demonstrating  several  features  of  the  method.  Both  problems 
have  also  been  solved  by  the  D-method  and  serve  as  a  basis 
for  limited  comparison  of  the  methods. 


I 

INTRODUCTION  TO  PROBLEMS  IN  OPTIMUM  CONTROL 

1.1  Formulation  of  a  problem 

We  desire  the  solution  to  the  set  of  ordinary  nonlinear 
differential  equations 

(1)  k±   =  fi[x1(t),x2(t), ,xn(t),p1(t),..,pm(t),t] 

i  =  1,2,... ,n;   to<  t  <  T 

which  is  optimal,  in  the  sense  of  maximizing  the  terminal 
value  of  a  performance  criterion 

s(xl> '^n*  t=T  ' 

while  simultaneously  satisfying  terminal  constraints 

(2)  \(Xi> i^&^tsT  =0  ;  k  =  1,2,. ..,N;   0  <  N  <  n. 

To  effect  these  conditions  we  have  at  our  disposal  m  con- 
trol or  driving  functions  p_(t),  s  =  l,2,....,m.   Their 
definition  for  W  t  ^  T  is  a  major  portion  of  the  problem. 


*  A  dot  (  )  over  a  variable  is  used  to  designate  the  deriv- 
ative with  respect  to  t. 

«t  If  an  integral  is  to  be  maximized  we  introduce  an  addi- 
tional state  variable,  x  +, ,  and  an  additional  differen- 
tial equation 


*n+l~  fn+l(xl>"->xn,t;) 


where  f  .,  is  the  integrand  of  the  integral 
n+j_ 


It  is  assumed  the  functions  f.   are  defined  over  an 
open  region  and  are  class  C" .   It  is  also  assumed  the  func- 
tions g  and  hk  are  independent  and  of  class  C1   in  the 
terminal  region.   The  variables  (x-,,x2, . . . ,  ,x  )   will  be 
termed  the  state  variables.   All  are  functions  of  the  in- 
dependent variable  t.  The  terminal  value,  T,  of  t  is 
either  specified  or  otherwise  determined  from  a  terminal 
condition.  The  values  (x-.  ,x2,  •  •  ,x  )  .   are  assumed  specified. 

It  will  be  convenient  to  represent  these  sets  of  varia- 
bles by  matrices  or  by  vectors.   We  define 


x  = 


x,(t) 

x2(t) 


xn(t) 


5  P  = 


P]_(t)" 

p2(t) 

• 

• 

p  (t) 

_ffl 

h  = 


l^Cxjt)' 
h2(x,t) 

hN(x,t) 


f  = 


fx(x,p,t) 
f2(x,p,t) 

« 

fn(x,p,t) 


No  distinction  will  generally  be  made  between  an  nxl 
(column)  or  a  l*n  (row)  matrix  and  a  vector.   If  a  quan- 
tity is  best  considered  a  vector,  as  in  a  scalar  product, 
it  will  be  designated,  for  example,  p. 
Equations  (1)  and  (2)  then  become 


(3) 
(4) 


x  =  f(x,p,t) 
h(x,t).  T  =0, 


We  desire  the  solution  of  system  (3)  which  maximizes 
the  performance  quantity  g(x,t)    subject  to  the  con- 

X> —  J. 

straints  (4) .  For  purposes  of  introduction  of  terminology 
and  some  basic  theory  a  brief  review  of  the  related  elements 
of  the  Calculus  of  Variations  is  first  presented. 

1.2  The  Classical  Calculus  of  Variations 

The  problem  formulated  in  section  (1.1)  is  typical  of 
those  considered  in  the  variational  calculus  -  those  of 


determining  maxima  or  minima  or,  in  general,  stationary 
values  (defined  on  page  8  )  of  functionals  such  as  g(x,t) , 
by  seeking  the  argument  functions  x(t)   and  p(t)   for 
which  the  functional  assumes  the  stationary  value  or  ex- 
tremum  in  question. 

An  outline  of  seme  phases  of  the  classical  method  of 
the  variational  calculus,  as  applied  to  a  greatly  simpli- 
fied optimum  problem,  will  serve  to  introduce  terminology 
and  to  present  essential  concepts  which  are  the  basis  for 
both  methods  of  solution  being  considered.   The  method  may 
be  termed  'indirect'  in  that  it  is  based  on  the  reduction 
of  the  variational  problem  to  differential  equations  the  so- 
lution of  which,  subject  to  the  boundary  conditions,  is  the 
solution  to  the  problem.   The  two  methods  considered  in 
this  study  may  be  termed  'direct'  in  that  they  consist  of 
construction  of  a  sequence  of  functions  that  converges  to 
the  desired  function. 

Following  Courant  and  Hilbert  [S],  we  consider  the 
simplest  problem  of  the  variational  calculus.   Find  the 
function  y(t)   which  is  such  that  y(t  )  =  A  and 
y(t1)  =  B,  Fig.  1,  and  the  functional 

,T 


F[y(t)]  =  f   f[t,y(t),y'(t)]  dt 


o 

is  maximized.   Any  function  y  =  y(t)   which  meets  the  end 
conditions,  is  of  class  C,  and  has  a  piecewise  continuous 


t    , 

o 


■■n(t) 


derivative  is     .en  ed     av.is- 
sible*      ■<-  t   there 

is  a   class  c£  admissible  argu- 
ments  and   that  among  them  there 
exists   the   solution   to  our 
problem,   the   function     y=y(t), 
which  maximizes   (or  minimizes, 
hereafter  we  will  refer  only 
to  maxima)    the  functional 
F  L  y  ( t )  j .      Su  p  p  o  s  e  we  hav  e 
found   the  desired   extremal  function     y  =  y(t).     Let  us  con- 
sider  also   an  unspecified  new  function     t)(t)      which  is  de- 
fined  for     t     n<   t, ,      possesses   the  necessary  derivatives 
o  x 

and   continuity,    and  is   such  that     Ti(tc)    =  rj ( t-^)    =  C,   but  is 
otherwise  arbitrary.     The  function 


Fig.  1  An  admissible 
variation  of  an  admissible 
curve 


y(t)  +  sri(t)  =  y(t)  +  6y(t), 


where  e  is  a  parameter,  represents  a  one-parameter  family 
of  admissible  functions  which  contains  the  solution  y(t) 
when  e=0«   The  quantity  6y(t)  =  e*n(t)  is  the  variation 
of  the  function  y  =  y(t).   For  sufficiently  small  magni- 
tudes of  €   the  varied  functions  y(t)   lie  in  an  arbi- 
trarily small  neighborhood  of  the  extremal,  y  =  y(t),  and 
the  integral 

F(y)  =  F[y(t)  +  eii(t)] 


may   be  regarded   as  a  function     i(c)       rhich  must  have   a  max: 
s  um   at     e  =  0      and   therefore  requires     i'(0)    =  0. 
Now 

y  =  y(t)   +  erKt) 
implies 

y'(t)    -  y'(t)    +   eti'(t) 
and  we  may  write 


$(e)    =    /      f[t,   y  +eri,   y'+eV]  dt 


o 

For  a  maximum  it  is  necessary  that 


(6)  i*(0)  =  /   nC  fy  -  |tfy,]  dt  =  ° 

O 

for  all  ri(t)  meeting  the  conditions  previously  imposed. 
The  fundamental  lemma  of  the  calculus  of  variations, 

[8]  page  185,  applied  to  (6)  requires 

or  equivalently, 

(7)  y"f  ,  ,  +  yT    ,  ,  +  f  ,   -  f   =  C. 

'      *  y '  y    "  y^1         y'x    y 

This  differed bj al  equation  was  first  discovered  by 
iuier  in  1744,  [9]  page  2>:,  and  is  commonly  referred  to  as 
the  Suler  equation.   A  solution  to  the  Buler  equation  is 
called  an  extremal.   Its  validity  is  a  necessary  condition 


for  the  existance  of  an  extremum;[8]  page  185.  Being  a 
differential  equation  of  the  second  order,  two  arbitrary- 
constants  must  be  determined  for  a  solution  which  meets  the 
boundary  conditions,  i.e.,  to  make  an  extremal  an  admissible 
extremal.   Every  solution  of  Euler's  equation  is  an  ex- 
tremal of  the  maximum  problem, [8],   Courant  describes  (6) 
the  variational  derivative  of  F  with  respect  to  t  and 
terms  its  role  here  as  analogous  to  that  of  the  gradient  in 
ordinary  maximum  problems. 

In  the  n+1  dimensional  case,  with  coordinates 

t,  y^t),  y2(t),  .  .  .  ,  yn(t) 


and   the  functional 

t 

•7lo 

the  Euler  equations  take  the  form 


tl 
I  =  /   f(t,  y1,  y2,  ...,  yn,  y£,  yg,..,  y^)  dt, 


(9)  «fy'  "  \   =  °5  ttV2  "  V  °5  -'M'  "  V  °' 

a  system  of  the  second  order  equal  in  number  to  the  number 
of  unknown  functions  yAt)    to  be  defined.   Any  n+1  space 

curve 

yl  =  yl(t)>  y2=  y2(t)>  »  yn=  yn(t) 

whi  ch  is  a  solution  to  the  system  (9)  is  an  extremal  and 
furnishes  a  stationary  value  of  the  functional  (8). 


8 


The  necessary  condition  of  Euler  is  usually  called  the 
first  necessary  condition.   There  are  also  necessary  con- 
ditions of  v/eiers trass,  Legendre  and  Jacobi  as  presented,  for 
example,  in  [l]  chapter  1.   These  are  not  discussed  here; 
the  fundamental  problem  confronting  an  applied  mathematician 
is  that  of  finding  a  likely  candidate  for  a  solution. 

1.3  Introduction  to  direct  numerical  methods 

The  problem  as  stated  in  section  (1.1)  is  similar  to 
that  resulting  in  system  (9).   The  state  variables  x(t)   of 
our  functional  g(x,t)  are  defined  by  a  set  of  n  first-order 
differential  equations  and  one  or  more  control  functions,  or 
parameters,   p(t);   t  is  the  independent  variable.  Solu- 
tions to  the  related  Euler  equations  will,  to  the  extent  of 
meeting  the  first  necessary  condition,  define  the  functions 
p(t)   such  that  the  resulting  path  x(t)  =  C  will  impart 
at  least  a  stationary  value  to  the  functional  g(x,t)m. 
The  proper  choice  of  the  arbitrary  constants  of  these  so- 
lutions will  result  in  admissible  extremals  if  they  are 
chosen  such  that  the  terminal  conditions  h(x,t)m  =  C  will 
also  be  met.   It  will  be  assumed  that  the  initial  conditions 

x(t)._.   are  specified  constants. 
z-zq 

Various  indirect  approaches  to  such  a  solution  exist. 
A  bibliography  to  indirect  method  studies  is  presented  in 
the  survey  paper  [10].  Analytical  solutions  of  the  Euler 
equations  involved  usually  require  idealizing  assumptions 


which  limit  their  applicability  in  practical  situations. 
Under  more  realistic  assumptions  a  numerical  attack  on 
these  equations  is  required.   Some  direct  methods  are  pre- 
sented in  [8]  page  174.   In  the  following  sections  the 
essentials  of  two  additional  practical  methods  are  out- 
lined.  As  with  all  methods  of  solving  Bolza's  problem, 
Lagrange  multipliers  are  employed  producing  equations 
equivalent  to  the  Euler  equations  in  terms  of  solutions 
to  an  adjoint  system  of  differential  equations.   The  con- 
cept was  developed  by  Bliss  [2]  in  connection  with 
differential  corrections  in  ballistics. 

The  D-method  uses  only  extremals.  Then,  by  an  iter- 
ative correction  procedure, obtains  an  admissible  extremal 
out  of  the  family  of  extremals. 

The  FL-method  produces  an  extremal  out  of  any  likely 
solution  to  the  original  equations;  part  of  the  problem  is 
in  a  subroutine  to  approach  an  admissible  solution.   The 
solutions  to  the  Euler  equations,  producing  the  extremal, 
are  approached  iteratively  by  a  gradient,  or  steepest- 
ascent  technique. 

Since  all  methods  of  solution  involve  solutions  to 
the  adjoint  equations  this  concept  will  be  outlined  be- 
fore the  details  of  the  two  methods  are  presented  seper- 
ately. 
1.4  The  adjoint  equations 

The  problem  as  stated  in  section  (L.  1)  was  to  determine 


10 


the  solution  to  the  system 

(3)  x  =  f(x,p,t) 

such  that  the  terminal  constraints 

(4)  h(x,t)T  =  0 

are   satisfied   and   the  functional  g(x,t)m  is  maximized. 

Of  primary  interest  is   the  relationship  between  vari- 
ations    6p(t)        in  the  control  variables     p(t),    at  any 
value     t,    and   the  resulting   terminal  variations     6xT     of 
the   state  variables.      We   also  need   the  relationships 
between     6p(t),    6huj,     and     6g™. 

Looking  at  (3)   in  component  formed),  the  variational 
relationships  for  the  state  variables  can  be  expressed   as 

*x<    =  Hi6xi   +  Hiftxo  +   •    •   HiSXn  +  ££i6Pi Hi5Pm 

i   9X,   1   8Xg  2       9x_  xi   9p,  rl  9p   ^m 

l  ■*  J»j  ^j  ••••)  n. 
In  matrix  notation  this  is 
(10)      6x  =  F  fix  +  G  6p 


11 


where 


F(t)    = 


dxi 


8£?  • 


9xi 


fi£n 

Ox*1 
n 


G(t)    = 


nxn 


8PJ   ' 


ap?    * 


Mi 
9Pm 


3  P. 


m 


J    nxm 


6p  = 


6p1(t) 


ePffi(t) 


6x  = 


6x, 


We  will  use  the  prime  symbol  to  denote  the  transpose  of  a 
matrix,   e0g«,  6x'  =  (fix-,,  6Xp,  .o...,  6X-) . 

Following  Bliss  [2]  the  system  of  linear  equations 
adjoint  to  (3)  is  defined  to  be 


(12) 


A=  -F'(t)A 


where 


t* 


A(t)  =  [X1(t),  X2(t) Xn<t)] 

is  a  matrix  of  n  Lagrange  multipliers,  and 

/  \  =  L  Ai )  ^J?>  •  o  o  «  •  •  ,  A_  J  • 

The  solutions  of  systems  (3)  and  (12)  are  related  by 


12 


dt[7\6x]    =  AG   6p; 


I  (-.  £}  7\ 


hence 

(13)  [ /Yfix   L.      =     /      (/\*G  6p)    dt. 

Since  we   are   assuming     x.      fixed     6x.    =  0  in   (13)    and  we  have 

o  o 

T 

(14)  [  Al6x   ]T    =    /t(A'G   6p)    dt 

o 

This  is   the  differential  formula  relating   the   terminal  changes 
in   the   state  variables   and   the  variations  in   the  control 
variables.      It  is   termed   the  Fundamental  Formula  by  Bliss   [2], 
It  is   also  known   as   the  Fundamental  Adjoint  Formula  and   as 
the  one-dimensional  form  of  Green's  Theorem  [11],    [12]. 

1.5  Applications  of  the  Fundamental  Adjoint  Formula 

Some  indication  of  the  usefulness  of  the  fundamental 
formula  is  displayed  in  the  following  observations. 

1)  If  A  is  chosen  as  the  specific  solution  to  the 
adjoint  system  (12)  such  that  at  t  =  T 


7\T   =  Vxxx(T)  =  (1,0,0,. ..,0), 


-**- 


■**-  7  x-,  refers  to  the  gradient  of  x-,  in   x  space. 


-x- 


We  assume  no  corners,  [9]  page  8,  exist, 
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let  us  call  this  solution  ^-(t),  (14)  then  becomes 

6xx(T)  =  J  A^   G  6p  dt  , 

to 

yielding  the  terminal  variation  of  x-,  due  to  variations 
6p(t).   Similarly  choosing  A^(t)   such  that 

A.(T)  =  Vx.(T)  =  (0,0, ...,1, ,0) 

i  =  1>2, ,n 

gives  m 

(15)  6x.  (T)  =  J.    A\    G  6p  dt  . 

1       zo     x 

2)   If  A  is  chosen  as  the  solution  to  the  adjoint 

system  such  that  A  (T)  =Vg(T),  (14)  implies 

g      x 

(16)  6g(T)  =  (Vg  *  6x)T  =  Jt  flg   G  6p  dt 

giving  the  first  variation  of  the  performance  function  in 
terms  of  the  parameter  variations  for  t  4-   t  ^  T. 


T 
(17)      611,(1)  =  J    A'j,  G  «P  dt  ,  i   =   1,  2,  ,  k. 


*o  i 


3)  For  the  constraint  functions  we  likewise  have 

.T 

4)  If  we  were  to  ignore  the  problem  of  meeting  con- 
straints and  look  for  a  solution  which  gives  a  stationary 
value  for  g(T),  we  would  require  the  vanishing  of  the  first 
variation, ( see  appendix  of  [4]) 

T 
(18)      6g(T)  =  J  /\'  G  6p  dt 
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which  requires     f\' G  =  C.      That  is 


^1,^2> 


»An)g 


Sl- 
ap: 


ap 


n 


8*1 

9  pi 


ef 


=  c 


which  represents  the  m  equations 


(19) 


A'.  •  §£  -  °-   i  -  i. 


,m. 


Equations  (19)  with  the  adjoint  differential  equations 
(12)  constitute  what  is  generally  termed  the  Euler-Lagrange 
equations  [4],  [13]  . 

In  equations  (18)  the  matrix  AG  is  clearly  an  in- 
fluence  function  giving  for  any  time  t  the  effect  that  a 
unit  impulse  fp?  at  that  time^wculd  have  on  the  perform- 
ance function  g(T).   In  that  regard  if  we  write  /\  G  as 
a  vector,  say  A(t),  (18)  then  can  be  written 


*«<*>   =  It 


A(t)  •  6p  dt 


and   for  maximum  change  in   the  performance  quantity     g(T) 


we  would   choose       fp        parallel  to    A(t)  for      t  <£   t  ^   T. 
It   should  be  noted   that     A(t)      and      fp(t)      are  m-dimensional 
vectors   ir    parameter   space   (p).      Courant,    [8]   page  223,   re- 
fers tc   vector     A(t)      as   the  gradient  of  the   functional 


g(  t) 


1IJ   x 


unction  space.   When  the  optimum  solution  or  path 
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in  function  space  is  obtained  this  gradient  will  vanish. 

In  the  case  of  our  general  problem  we  have,  in  addition 
to  the  requirement  of  maximizing  the  functional   g(x,t)m,  the 
necessity  for  meeting  terminal  constraints  h(x,t)  =  0.  We 
will  choose  the  vector  p(t)   in  parameter  space  such  that 
the  path  meeting  the  constraints  yields  the  largest  possible 
value  for  g(T).   We  now  consider  two  methods  of  solving  the 
problem. 

We  have  called  g(x,t)m  a  functional.  It  is  a  function 
which  becomes  a  functional  since  (x,t)  represents  the  endpoint 
of  a  curve^  satisfying  the  given  differential  equations  with 
initial  point  x(t  )   given.  The  function  g  has  become  a 
functional  in  the  sense  that  the  point  (x,t)  depends  upon 
the  curve. 
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II 

TV/C  NUMERICAL  METHODS 

2.1        The  differential  method 

Recalling   the  fundamental  differential   formula 

(14)  [X.x]T     =/tT/\'G6pdt    , 

o 

the  D-method  is  based  on  the  Euler  equation  which  states 
that  if  the  path  is  to  furnish  a  maximum  to  some  functional 
g(x,t)m,  it  is  necessary  that  the  coefficient  Ag  of  6p 
in  (14)  vanish  for  some  solution  A  of  the  adjoint  system 
for  all  t  ^  t  ^  T.  This  requirement,  Ag  =  0,  expanded 
in  terms  of  (11),  implies  that  for  all  tj£   t  ^  T  there 
exists  a  A(t)   and  x(t)   such  that 

Att)  .  J[  =  o 


Act)  .   |f  =o 

(20)  Sp2 


Act)  .  |f  =o 

8Pm 

The  problem  becomes:  determine  A(t)   and  hence  the  control 
functions  p(t)   such  that  (20),  (3),  (12)  and  (4)  are  satisfied, 

Assuming  first  that  T  is  specified  we  proceed  to  out- 
line the  method  of  determining  the  A  matrix.  On  page  14 
a  set  of  n  solutions  of  the  adjoint  system  (12)  was 
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developed  such  that 

A±(T)  =  VXi(T),   i  =  1, ,n. 

The  linear  combination 

(21)      A(t)  =  c1/\1   +  c2/\2  +  +  cn^n 

of  these  n  solutions  is  also  a  solution  to  (12). 

As  a  first  approximation  to  the  admissible  extremal 
solution  desired  we  estimate  a  set  of  constants  c.   for 
(21)  and  determine  the  associated  extremal  by  employing 
this  A(t)   in  system  (20).   This  will  in  general  not  be 
an  admissible  extremal  as  it  will  not  meet  the  terminal 
constraints  (4) . 

On  the  basis  of  the  error  in  meeting  the  terminal  con- 
ditions j-  corrections  dc^  for  the  constants  in  (21)  are 
determined.  The  resulting  improved  A(t),  when  employed 
with  (3)  to  obtain  a  new  solution,  will,  if  the  method  con- 
verges, give  a  more  nearly  admissible  extremal.   The  process 
is  repeated  until  some  convergence  criterion  is  satisfied. 

If  the  terminal  T  is  not  specified  it  too  will  be  es- 
timated for  the  first  solution  and  a  correction  dT  will  be 
produced  with  the   dc.  at  the  end  of  each  iteration. 

A  more  detailed  analysis  of  the  basis  and  the  calcu- 
lation procedure  is  presented  in  appendix  (A). 
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2.2  The  method  of  the  fundamental  lemma 

This  method  is  concerned  directly  with  the  construction 
of  the  set  of  m  control  variables  p(t),  for  t  ^  t  ^  T, 
which  produce  the  solution  to  (3)  that  satisfies  the  terminal 
conditions  (4)  and  maximizes  the  performance  quantity 
g(x,t)j.  Essentially  the  procedure  is  to  produce  any  reason- 
able trial  solution  by  the  choice  of  a  nominal  control  pro- 
gram  p(t).   The  resulting  solution  will  in  general  be 
neither  admissible  nor  an  extremal.  That  is,  the  con- 
straints (4)  are  not  satisfied  and  the  Euler  equations 
(19),  a  necessary  condition,  are  not  satisfied,,   A  better 
solution  can  therefore  be  constructed,,   A  variation  6p(t) 
is  then  produced  such  that  the  Euler  equations  are  more 
nearly  satisfied  and  the  solution's  deviation  from  that 
of  an  admissible  extremal  is  reduced ,  This  variation  ap- 
plied  to  the  nominal  p(t)   produces  a  better  control  pro- 
gram ptt)  =  ptt)  ■*  6p(t)e  The  procedure  might  be  de- 
scribed as  iteratively  stepping  a  control  vector  p(t), 
defined  in  parameter  space  for  t  <  t  >$  T,  in  the  direc- 
tion of  steepest  ascent  toward  meeting  constraints  and 
toward  meeting  the  first  necessary  conditions  of  an  ex- 
tremal. This  direction  is  that  of  -V  h  and  V  g  in 

P        P 

parameter  space,  where  we  use  the  subscript  (p)  to  denote 
parameter  space.   As  an  admissible  solution  approaches  an 
extremal  V  g  will  tend  to  zero  [6]c   The  direction  of 
6p,  for  each  point  of  a  nominal  solution,  will  be  determined 
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by  the  direction  cf  the  gradients  at  that  point.   The 
magnitude  of  6p  at  each  point  will  be  best  chosen  pro- 
portional to  the  magnitude  of  the  gradients  at  that  point 
in  parameter  space. 

In  appendix  B  the  desired  parameter  correction  program 
p(t)   is  derived  and  defined  for  t  <  t  <  T  as 

( 22)  J 

,      /(ds)^-  dh'lj-  dh  , 

6P(t)  =iG-(Ag0-Ahn  ij£  ihg)/- f   +  G'AhQ  1-1  dh  . 

j  igg"  lhg  lhh   xhg 

The  plus  sign  is  used  if  g(x,t)T  is  being  maximized. 

One  phase  of  the  calculation  procedure  as  outlined  by 
Bryson  [6],  and  summarized  in  appendix  B,  was  found  to  be 

somewhat  indefinite.   It  is  that  of  determining  the  arbi- 

p 

trary  constant  (ds)   which  governs  the  magnitude  of  the 

step  size  6p(t)   by  limiting  the  integral  of  these  mag- 
nitudes for  t<  t  <  t  by  the  relationship  (B-7)  .   It 

p 

seems  apparent  that  the  larger  the  magnitude  (ds)  ,  within 

the  limit  which  allows  the  linearization  (10),  the  more 

rapidly  our  solution  process  will  converge  to  the  admissible 

2 

extremal.   Accordingly  the  first  restraint  on  (ds)   will 

2 

be  imposing  an  upper  limit  (ds  )   as  estimated  from  the 

nature  of  the  problem  in  conformance  with  (10). 

Further  it  is  obvious  that  repeated  application  of 
large  magnitude  variations  ^p(t)   will  not  allow  fine  ad- 
justments toward  the  admissible  extremal,  as  the  solution 


o, 


0 


is  approached,  unless  the  magnitude  (ds)^  is  further  ma- 
nipulated.  We  therefore  desire  an  automatic  scheme  which 
will  rapidly  step  the  nominal  solution  toward  the  ultimate 
solution  but  will  adjust  itself  to  finer  increments  in  the 
vicinity  of  the  admissible  extremal.   The  following  dis- 
cussion of  Bryson's  procedure  is  to  develop  a  basis  for  the 
design  of  such  an  automatic  convergence  scheme. 

The  expression  (22)  defining  the  vector  6p  consists 
of  two  components  which  we  define  6p,   and  {p  ,  We 
designate 

(23)       5ph=  G'A^I^dh 

since  it  is  essentially  a  vector  in  p  space  whose  magni- 
tude is  proportional  to  the  magnitude  of  the  error  in  meet- 
ing the  constraints  resulting  from  the  nominal  solution. 
6p,   has  the  direction,  essentially,  of  -V  h.   Accordingly 
6 p.   is  a  component  of  ftp  which  produces  a  more  nearly 
admissible  solution.   The  desired  correction  toward  an 
admissible  solution  may  be  large  and  may  require  larger 
values  of  lap  I  than  the  upper  limit,  previously  imposed, 
will  allow.   We  give  construction  of  an  admissible  so- 
lution first  priority,  however,  and  up  to  the  limit  (ds  ) 
demand  that  gp  be  designed  to  meet  the  constraints. 

The  first  term  of  (22),  involving  the  radical,  is 
essentially  a  vector  in  p-space  whose  direction  is  normal 

to  V  h  and,  as  nearly  as  meeting  constraints  will  allow,  is 

P 

parallel  to  V  g,  It  specifies  a  component  6p  which  causes 
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progression  toward  maximization  of  g(x,t)T  without  dis- 
turbing the  progress  toward  an  admissible  solution  made  by 
6 p..   The  term 


'(ds)2-  dh'  Ij£  dh 


2 

is  what  remains  of  magnitude  (ds)   for  use  as  magnitude  for 

a  vector  normal  to  opy.*   If  6p,   demands  all  of  (ds) 
this  radical  is  zero  and  6p  =  6 p..  Progress  will  be  made 

toward  an  admissible  solution  and  the  functional  g(x,t)m 

will  increase  or  decrease  according  to  whether  or  not  ~^Pu 

has  a  component  in  or  opposite  to  the  direction  of  ^Dg« 

Accordingly,  the  scheme  for  the  choice  of  the  magni- 

(ds)2  will  1 
limit,  a  magnitude 


2 
tude  (ds)   will  be  to  first  provide,'  within  the  upper 


I  G'Ahn  *Sh  dh  I 

for  attaining  an  admissible  solution,  and  second,  provide 
for  a  component  toward  the  extremal,  adjusted  according  to 
the  estimated  deviation  from  an  extremal .   We  define  this 
latter  magnitude  (A) ,  the  magnitude  of  6p  .  We  then  re- 
quire 

(24)      (dso)2  >y     (ds)2=  A2  4    |G»  AhQ  Ij£  dh  | 

The  initial  choice  and  automatic  adjustment  of  the 
magnitude  (A)  is  another  problem.   We  want  it  as  large  as 
possible  initially.   This  will  be  controlled  by  (24). 
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jv'e  want  smaller  magnitudes  as  the  ultimate  solution  is 
approached,  am torn a tic ally  adjusted  to  the  need.   The  con- 
vergence criterion  suggested  by  Bryson  [6]  v/ill  now  be  em- 
ployed.  He  develops  the  fact  that  on  an  admissible  solu- 
tion the  denominator  of  the  radical  of  (22)  is  the  direc- 
tional derivative  or  gradient 


(25)      df  "Jlgg  '   Jhg  Jhh  Jhg 

in  p-space.   This  must  decrease  as  the  extremal  is  approached. 
This  quantity  is,  however,  in  general  greater  than  zero.  It 
is  reasoned  that  if  a  relatively  large  (A)   is  repeatedly 
applied  to  an  admissible  solution  it  will  rapidly  decrease 
the  quantity  dg/ds  of  (25)  and  then  will  ultimately  cause  it  to 
pass  through  a  minimum,  greater  than  or  equal  to  zero.   If 
at  this  point  the  magnitude  (A)   is  not  reduced  subsequent 
solutions  will  be  over-corrected  toward  an  extremal;  trajec- 
tories which  oscillate  about  the  solution  to  the  problem. 
We  use  this  characteristic  to  trigger  a  reduction  of  (A) 
to  some  fraction,  say  0.2  A,   and  step  from  there  toward 
the  extremal  in  smaller  steps 0  The  minimum  value  for  dg/ds 
attained  this  time  will  be  smaller  and  the  same  behavior  will 
signal  the  need  for  further  reduction  in  (A) .  The  process  is 
repeated  until  some  lower  limit  on  (A)  or  dg/ds 9  a  conver- 
gence criterion,  is  attained0 

A  flow  chart  of  a  convergence  scheme  is  presented  as 
appendix  C0 
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Ill 

APPLICATIONS  AND  COMPARISON 

3.1  A  basis  for  comparison 

The  performance  of  both  methods  was  investigated  for 
two  ship  routing  problems.   Of  primary  interest  were: 

(a)  the  problem  of  choice  of  starting  para- 
meters that  result  in  convergence  to  a 
solution. 

(b)  the  domain  of  these  parameters  that  pro- 
duces a  solution. 

(c)  the  range  of  problem  variations  that  can 
be  solved  with  an  acceptable  set  of  these 
parameters. 

(d)  the  rate  of  convergence  to  a  solution. 
For  each  application  the  performance  criterion  was 

defined  as  the  functional 

_  (x  -  t/100)2+  4v2 
(26)      g  =ftT   [  c  +  e  2        ]  dt 

which  could  be  interpreted,  for  example,  as  a  probability 
of  detection  or  as  a  measure  of  accumulated  radioactive 
fallout  in  an  area  where  the  intensity  is  defined  by  the 
time  and  position  function  g(x,t)e   This  is  essentially 
the  problem  considered  by  Faulkner  [4]  and  Cook  [14]. 
Their  programs  for  the  D -method  were  employed  as  a  source 
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of  comparison  data  for  that  method. 

More  specifically  the  problems  considered  are: 
Problem  I:  Determine  the  route  from  a  point  within  an 

area  where  £(x,t)   is  defined,  to  any  point 
at  which  g(x,t)   is  a  prescribed  tolerable 
maximum,  such  that  the  accumulated  fallout 
g(x,t)T  is  minimized. 
Problem  II:  Determine  the  route  between  two  prescribed 
points  in  the  area  where  g(x,t)   is  defined 
such  that  g(x,t)T  is  minimized. 
For  the  D-method  the  calculation  procedure  of  page  18 
for  these  problems  is  reduced  to  the  specification  of  a 
pair  of  constants  defined  by  Cook  [14]  as  (X  ,u  )   and  are 
the  two  components  of  the  initial  adjoint  vector,  which  are 
unknown  at  the  outset  of  computation. 

For  the  FL-method  the  starting  parameters  are,  as  for 
any  problem,  a  nominal  p*(t)   and  a  magnitude  (ds  ).   For 
these  problems  p*(t)   is  a  bearing  angle  or  heading  and 
could  be  any  constant  direction  p*(t)  =  k  radians,  rough- 
ly away  from  the  point  of  maximum  g(x,t),  'ground  zero1. 

From  (B-7)  an  acceptable  value  for  (ds  )  is  determined  by 

2        2 

(ds  )   =  (6p)  T,  where  T  is  approximately  the  total  time 

required  and  6p  is  the  change  in  bearing,  at  any  point, 
that  would  reasonably  meet  the  linearity  requirements  of 
(10).   As  an  example;  estimating  T  as  100  with  an  allow- 
able 6p  of  0.5  radians  would  give: 
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(dso)*  =   (C.5)"(10C)  =  (5)*, 


3.2  Problem  I 

We  consider  the  specific  case  of  problem  I  with  start- 
ing point  (*o,yo)   specified  as  (0.3,0.1).   For  the  FL- 

methcd  the  choice  dS  =5  is  reasonable  and  a  logical 

o 

choice  for  p*(t)   is  p*(t)  =  tt/2  for  C  <  t  ^  T  or  very 
roughly  'away'  from  g(max)   at  (0,0).   For  the  D -method  a 
reasonable  choice  (with  the  advantage  of  hind-sight)  for 
(X  ,u  )  is  (0,-100)  which  is  a  vector  roughly  in  the  direc- 
tion of  the  negative  gradient  of  g(x,t)   with  a  magnitude 
approximating  the  final  time  T. 

In  an  effort  to  experimentally  determine  the  intervals 
of  convergence  an  interval  of  values  containing  these  esti- 
mated starting  parameters  was  explored  for  both  methods. 
The  results  are  summarised  in  Table  I  and  Figures  2  and  3. 
For  the  FL-method  the  process  converged  for 
0.1  <  p*(t)  ^  3.2  radians.   The  starting  parameter  is 
therefore  any  bearing  angle  roughly  'away*  from  (0,0)  within 
about  a  180°  interval.   In  Fig.  2  path  A  depicts  the  op- 
timum route.   For  values  in  the  vicinity  of  3tt/2  the  local 
minimum  at  the  opposite  side  of  the  cloud  (path  B  in  the 
sketch)  is  developed.   If  headed  in  the  direction  p*(t)=  v y 
however,  the  local  maximum  will  not  be  produced  but  the 
nearest  minimum  (path  A)  will  be  generated,  while  p*(t)=  3.6 
results  in  convergence  to  pathB. 
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Table  I   Starting  parameters  for  problem  I. 


D-Method 

FL-Method 

>o 

no 

Convergence 

P~(t) 

Convergence 

90 

-10 

yes 

0.1 

yes 

80 

-20 

yes 

0.2 

yes 

60 

-40 

yes 

0.4 

yes 

40 

-60 

yes 

0.8 

yes 

20 

-80 

yes 

1.2 

yes 

0 

-100 

yes 

1.6 

yes 

-20 

-80 

no 

2.0 

yes 

-40 

-60 

no 

2.4 

yes 

-60 

-40 

no 

2.6 

yes 

-80 

-20 

no 

3.0 

yes 

-90 

-10 

no 

3.2 

yes 

0 

-50 

no 

3.6 

(path  B) 

0 

-200 

yes 

60 

40 

(path  B) 

80 

20 

ii 

0 

100 

it 

-60 

40 

no 

-80 

20 

no 
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converge  to  path  A 
P  =  3.2^ 


p  =  3.6  *z  " 
converge  to  path  B 


Figure  2   FL  Convergence  Interval  for  p*(t)=k 


The  number  of  iterations  required  for  a  'solution'  de- 
pends on  the  accuracy  demanded.   For  values  within  the  inter- 
val of  convergence  the  poorest  choice  for  p*(t)  requires 
about  ten  percent  more  iterations  than  required  of  the 
best  constant  choice.  No  great  penalty  is  associated  with 
rough  estimation  of  the  starting  parameter,  for  a  close 
approximation  to  an  admissible  route  is  produced  rapidly. 
Computer  time  is  of  the  order  of  45  seconds  for  40  itera- 
tions.  The  first  15  iterations  produced  a  very  accurate 
solution. 

For  the  D-method  starting  parameters  were  explored  in 
the  intervals  -100^  xo  ^  100,   -100  <  Mo  <  100°  For  Pur- 
poses of  probing  the  influence  of  magnitude  (Xo,uo)  =(0,-50) 
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and  (0,-SOC)  were  also  tried.   Figure  3  and  Table  I  sum- 
marize the  results.   In  each  case  for  which  a  solution 
existed  convergence  was  relatively  rapid  requiring  about 
10  iterations  and  comouter  time  in  the  order  of  10  seconds 


no  solution 


converge  to  path  B 


(90,-10) 


converge  to 
path  A 


Figure  3  D-method  Convergence  Interval  for  (X  ,u  ) 


for  the  particular  convergence  criterion  employed. 

Continuing  with  problem  I  we  explore  the  domain  of 
starting  points  (x  ,y  )  for  which  a  given  set  of  starting 
parameters  produces  a  solution.   For  practical  purposes  it 
would  be  desirable  that  a  given  set  of  parameters  func- 
tion for  a  large  domain  of  starting  points  under  the  radio- 
active cloud.   Accordingly  solutions  were  attempted  for  a 
number  of  starting  points  in  an  interval  bounded  roughly  by 
the  left  and  right  extremeties  of  the  intolerable  fallout 
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area.   The  input  parameters  were  maintained  at  (X  ,u  )=(0,-100) 
for  the  D-method,  and  p*(t)=1.0  for  the  FL-method. 

Figure  4  and  table  II  portray  the  results  for  eight  start- 
ing points.  Both  methods  produced  solutions  for  each  rune  In 
figure  4  the  terminal  points  are  associated  with  the  corre- 
sponding initial  point  by  a  dashed  line  simulating  the  route. 
As  table  II  indicates,  the  coordinates  of  the  terminal  points 
agree  for  the  two  methods  to  the  third  or  fourth  significant 
figure.   Values  for  final  time  (T)  and  the  performance  g  like- 
wise agreed  to  at  least  four  figures.  Extended  results  for  them 
are  given  only  for  run  D  , 

In  regards  to  the  FL-method,  table  III  is  a  printout  of 
the  terminal  values  of  the  iterations  toward  a  solution  for 
run  D  of  figure  4  and  table  II.  It  will  serve  to  display 
several  features  of  the  FL-method  as  applied  to  a  relatively 
extreme  choice  of  the  starting  parameter  p*(t).  Specifically 
the  starting  point  is  (0.3,0.1)  and  p*(t)   was,  for  probing 
purposes,  chosen  as  p*(t)=0.1.   This  is  an  obviously  unwise 
choice  in  as  much  as  the  cloud  and  the  ship  have  comparable 
speeds  and,  for  this  case,  are  almost  parallel.  Consequently 
the  terminal  time  (T)   for  the  nominal  run  is  about  ten 
times  the  value  on  the  extremal.   A  nominal  solution  is  pro- 
duced, however,  and  from  it  relatively  rapid  convergence  to 
an  extremal  is  portrayed. 

The  parameter  (ds )  was  chosen  as  (6.0).   The  choice  of 
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Fig.   4     Allowable   Starting  Points   for   (>c,m0)    and  p   (t)    fixed 


Table  II   Allowable   Starting  Points   (xo,yQ)    for  Problem  I 
(>o,u0)    =   (0,-100);        p*(t)    =  1.0 


Start  point 
(xo>yo} 

X(T) 

Y(T) 

FL 

D 

FL 

D 

A  (1.2, .1) 

1.4092 

i 

1.4103 

1.0511 

1,0510 

B  (1.0, .1) 

1.1379 

1.1388 

lc0702 

i 

1.0702 

C  (0.7,0.1) 

0.7266 

0.7223 

1.0661 

1.0662 

D  (0.3,0.1) 

0 . 1775 

0.1778 

1.0074 

1.0C75 

E  (-0.7,0.1) 

-1.1066 

-1.1071 

0.6128 

0.6125 

F  (-1.0,0.1) 

-1.4175 

-1.4179 

0.4406 

0.4401 

G  (-1.2,0.1) 

-1.5879 

-1.5882 

0.3325 

0.3320 

H  (-1.4,0.1)  | 

-1.7292 

-1.7293 

0 . 2454 

0 . 2452 

Typical  values   for   final   time  T   and   performance   g(T) 
for   point  D:      T^T    =  91.6399,    Tn   =  91.6417 


'FL 


g(T)F]   ■=  49.0500,    G(T)D   =  49.0501 
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Table  III 

Problem 

I,  FL-Method,  Run  D 

X^  =  0,3000    Y^  = 
o              0 

0 . 1000 

S  =  6,00 

pX(t)  = 

0.10 

T 

X(T) 

Y(T) 

g(x,t)T 

D 

ds 

967,1896 

9.9236 

1.0656 

487.8577 

17661.3315 

6.0000 

256 . 7625 

2.6762 

1.0716 

127,6031 

2228 . 5757 

6.0000 

124.8307 

1.0787 

1.0696 

62.3428 

487.7631 

6.0000 

95 e 2305 

.3507 

1.0300 

49.8881 

79 . 2727 

6.0000 

98.7530 

.1648 

.9910 

53.9286 

275.8876 

1.8000 

94.0974 

.1323 

.9939 

50.3492 

132.3472 

1.8000 

91.8374 

.1091 

.9938 

49 . 1760 

29.3284 

1.8000 

92.2926 

.2244 

1.0146 

49,3411 

64.2659 

.5400 

91.8649 

.2062 

1.0121 

49.1077 

23.8514 

.5400 

91.6341 

.1808 

1.0080 

49,0658 

15.8866 

.5400 

91.7693 

.1826 

1.0081 

49 . 1008 

28.9107 

.1620 

91.6957 

.1801 

1.0078 

49.0645 

15,1434 

.1620 

91.6437 

.1775 

1,0074 

49,0503 

1.5992 

.1620 

91.6199 

.1759 

1.0072 

49.0569 

11.0875 

.0486 

91.6262 

.1765 

1.0073 

49.0528 

7.1449 

.0486 

91.6340 

.1771 

1.0074 

49.0505 

3.2402 

.0486 

91.6426 

.1776 

1.0074 

49,0501 

.7613 

.0486 

91.6293 

.1766 

1.0073 

49,0503 

3.1283 

.0146 

91.6331 

.1769 

1.0073 

49.0501 

1,9116 

.0146 

91.6372 

.1773 

1.0074 

49.0500 

.7058 

.0146 

91.6416 

.1776 

1.0074 

49.0501 

.4999 

.0146 

91.6376 

.1773 

1,0074 

49,0500 

.7230 

.0044 

91.6388 

.1774 

1.0074 

49.0500 

.3542 

.0044 

91.6399 

.1775 

1.0074 

49.0500 

.0208 

.0044 

91.6395 

.1774 

1.0074 

49,0500 

.3385 

.0013 

91.6396 

.1774 

1.0074 

49,0500 

c  £j£-0^£ 

.0013 

91.6397 

.1774 

1.0074 

49.0500 

.1132 

.0013 

91.6399 

.1774 

1.0074 

49.0500 

.0149 

.0013 

91,6398 

.1775 

1.0074 

49.0500 

.0327 

.0004 

91.6399 

.1775 

1,0074 

49,0500 

.0116 

.0004 

91.6399 

.1775 

1.0074 

49,0500 

,0062 

.0004 

91.6398 

.1775 

1.0074 

49.0500 

.0189 

.0001 

91.6399 

.1775 

1,0074 

49.0500 

.0106 

.0001 

91.6399 

.1775 

1.0074 

49.0500 

.0031 

.0001 

91 . 6399 

.1775 

1.0074 

49.0500 

,0040 

.0000 

91.6399 
Computer 

.1775 

times  44 

1.0074 
seconds 

49,0500 

,0000 

.0000 
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any  number  0  <  dsQ<  20  would  also  produce  solutions .  Actually 
values  nearer  20  would  hasten  the  convergence  while  values 
smaller  than  six  would  require  a  greater  number  of  iterations. 
Table  III  also  portrays  the  operation  of  the  convergence 
scheme  discussed  on  page  23 „   Since  in  this  case  there  are  no 
constraints,  A  in  (24)  is  (ds)0   As  table  III  indicates,  dS 
stays  at  its  initially  chosen  value  of  (660)  for  four  itera- 
tions o  By  this  time  the  extremal  has  been  closely  approxi- 
mated not  withstanding  the  difficult  nominal  route  chosen.   For 
the  fifth  iteration  dg/ds  has  crossed  through  its  minimum  and 
the  reduction  of  (ds)  to  0.3(dS)  has  been  switched.   The  pro- 
cess then  continued  until  an  extremal  was  crossed  again, caus- 
ing a  further  reduction  in  step  size0  In  this  case  the  pro- 
cess  was  halted  when  the  condition  dg/ds  <  10  was  attained. 

3.3  Problem  II 

The  two  methods  were  less  directly  compared  for  problem  II. 
Experimental  probing  was  directed  primarily  at  weaknesses  of 
the  convergence  scheme  devised  for  the  FL-method.   The  first 
specific  problem  of  type  II  was?  Determine  the  optimum  route 
from  (xo,yo)=(0o03,0.1)   to  (xf ,yf )=(lo0,lo0)   with  c  =  0 
in  (26).  An  admissible  solution,  simulated  by  route  (A)  in 
figure  5,  was  closely  approximated  in  the  third  iteration . 
Convergence  to  an  admissible  extremal,  however,  was  not  re- 
alized.  Longer  and  longer  admissible  routes,  (B),  were  devel- 
oped, each  better  than  the  previous  from  the  standpoint  of 
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decreasing  the  performance  criterion,  g  .  It  then  became  appar- 
ent that  the  solution  to  this  problem  does  not  exist.  The  quan- 
tity g  clearly  would  be  minimized  by  escaping  from  under  the 
radioactive  cloud  to  some  very  distant  point,  returning  after 
the  cloud  has  moved  away.   The  FL-method  developed  a  sequence 
of  solutions  which  approached  this  solution.  The  D -method  for 
this  case  returned  an' error'  printout  in  as  much  as  one  of  the 
input  parameters  is  a  relatively  close  estimate  of  the  final 
time  T.  This  value  being  undefined  the  process  stopped  due 
to  too  large  an  error  in  the  estimate  of  the  magnitude  of  this 
parameter. 

The  FL-method  was  then  applied  to  the  same  problem  with 
c  =  0.1  in  (26).  An  admissible  extremal  resulted  in  a  reason- 
able number  of  iterations. 

The  importance  of  the  careful  choice  of  the  stopping  con- 
dition 0(x,t)m=  0,  for  the  FL-method,  was  demonstrated  by 
problem  II  for  the  route  (0.03,0.1)  to  (0. 0, J. 0) ,.  figure  <5, 


'  ~>x  B 

{f     *"tt  1Y 


k 


(.03,.l) 


g(t)=k 

\ 


/ 


i(t)=k 


(.03,.l) 


Figure  5  Problem  Ila 


Figure  6  Problem  lib, He 
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with  fi(x,t),p=  0  chosen  as  (x  -  xF)m  =  0,   The  resulting 
situation  is  that  the  stopping  constraint  is  nearly  parallel 
to  the  route  making  T,  and  thus  x(T)  and  g(x,t)T  particularly- 
sensitive  to  variations  in  the  control  variable,,  Consequently 
difficulty  arose  in  that  the  constraint  (x  -  xF)m  =  0  was 
difficult  to  meet  with  the  accuracy  desired.  The  choice  of  the 
line  (y  -  yF)T  =  0,  or  the  circle  through  (xF,yF)  with 
center  (x  >y  ),  would  alleviate  the  situation,,   In  general 
n(x,t)^  =  0  and  g(x,t)T  =  0  should  be  as  nearly  parallel 
as  possible. 

The  difficulty  arising  from  too  large  a  value  for  A  was 
demonstrated  for  the  FL-method  by  the  solution  of  problem  II 
for  the  route   (0.03,0.1)  to  (-1,0.1),  simulated  by  curve  (D) 
of  figure  6.   The  initial  choices  of  magnitudes  ten  and  five 
for  A  produced  very  slow  convergence  toward  an  admissible 
route.   A  =  0.1,  however,  converged  reasonably  well.  This 
can  apparently  be  attributed  to  errors  due  to  linearization 
on  a  path  wi th  relatively  large  curvature. 

The  D-method  produced  rapid  convergence  to  solutions  for 
each  of  the  last  two  cases. 

3.4   Conclusions 

On  the  basis  of  this  study  the  following  observations 
regarding  the  points  of  comparison  (a)   through  (d)  of 
section  3.1  seem  justified. 
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(a)  The  choice  of  acceptable  starting  parameters  is  di- 
stinctly easier  for  the  FL-method0   The  matrix  p(t)   is  the 
set  of  actual  control  parameters  for  the  problem  and  as  such 
it  is  usually  a  set  of  well  understood  physical  quantities,, 
These  are  less  abstract  than  the  adjoint  vector  or  the  con- 
stants C.  of  the  D-methodo   General  knowledge  about  the 
physical  concepts  involved  usually  permits  estimation  of  a 
starting  control  program  which  relatively  closely  approximates 
that  of  the  optimum  solution. 

(b)  The  interval  of  convergence,  i6e0,  the  domain  of 
starting  parameters  from  which  a  solution  can  be  produced,  is 
significantly  larger  for  the  FL-method9   Any  reasonable  esti- 
mate for  the  control  parameter  will  start  the  solution .  Table 
I  and  figure  3  indicate  there  are  seemingly  good  estimates 
for  (*0>u 0) >  such  as  (-20,-80)   which  do  not  produce  solu- 
tions for  the  D-method0 

(c)  Once  a  workable  set  of  starting  parameters  is  found, 
however,  both  methods  produce  solutions  from  it  for  a  wide 
range  of  variations  of  the  problem . 

(d)  If  the  D-method  converges  it  consistently  results  in 
the  more  rapid  convergence.   For  similar  convergence  criteria 
the  computation  time  ratio  is  about  one  to  two,  and  the  number 
of  iterations  required,  a  ratio  of  about  one  to  four,,   The 
criteria  employed  were  different  for  the  two  methods,  however, 
and  it  appears  from  table  III  that  the  requirement  dg/ds  <  10 
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was  too  tight  for  practical  purposes.   The  precision  is  avail- 
able, however . 

The  FL-method  took  only  about  ten  percent  more  iterations 
to  converge  from  a  relatively  wild  estimate  for  p*(t)   than  it 
took  for  a  very  close  estimate „ 

The  following  comments  regarding  the  two  methods  are 
based  on  experience  in  implementing  the  two  methods  and  on 
analyses  of  the  papers  by  Faulkner  [4]  and  Bryson  and  Denham 
[6]. 

The  requirement  of  storing  an  entire  control  program  and 
a  complete  set  of  functions  A'  G  for  any  one  iteration  is  a 
disadvantage  for  the  FL-method  .   This  could  result  in  memory 
storage  problems  for  large  problems  or  small  computers <,  In 
contrast  the  D=method  holds  only  the  parameters  C.  and  T  or 
their  equivalent 0 

As  indicated  in  Appendix  B  the:. FL-method  has  been  develop- 
ed to  the  point  of  providing  for  variations  at  t=t  and  for 
the  insertion  of  weighting  functions  by  the  matrix  W(t)0  The 
power  of  these  features  has  not  been  investigated  in  this 
paper. 

The  problem  of  discontinuous  control  (corners)  or  un- 
bounded control  have  also  been  ignored 0  The  D-method  has  been 
developed  to  handle  them?  [4],  while  the  application  of  the 
FL-method  to  such  problems  has  not  been  undertaken,  nor  is  it 
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clear  how  this  could  be  dene. 

With  the  D-method  there  is  nc  question  of  how  closely  we 
have  approximated  an  extremal, with  the  mathematical  routine, 
since  only  extremals  are  used.   The  only  limitation  is  the 
accuracy  of  the  integration  method 8   Also,  since  only  extremals 
are  used  various  identities  associated  with  an  extremal  can 
be  employed.   For  example,  if  the  independent  variable  does 
not  enter  explicitly,  the  Hamiltonian  is  constant. 

The  integral  relations  for  the  corrections  in  the  FL- 
method  seem  to  be  more  stable  from  the  computational  standpoint, 
while,  if  the  strong  Legendre  condition  is  not  satisfied  the 
routine  must  be  modified  in  the  D-method . 

Only  forward  integration  (or  only  backward  integration) 
is  used  in  the  D-method .  In  the  FL-method  the  basic  equations 
are  integrated  forv/ard  and  the  corrections  are  obtained  by 
backward  integration  of  the  adjoint  systemc  This  necessitates 
either  storing  the  complete  trajectory  just  produced  or  re- 
constructing it  from  terminal  values  by  backwards  integration 
of  the  basic  differential  equations  of  the  system . 
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APPENDIX   (A)   THE  DIFFERENTIAL  METHOD 

A.l  Outline  of  the  analysis 

The  following  relations  have  been  defined: 


(3) 


x  =  f(x,p,t)   the  state  variable  equations 


(4)       h(x,t)m       N  terminal  constraints. 


(12) 


a   =  -F'A 


the   adjoint   system, 


(14)  [A'6x]T  =    I     A1    G  6p  dt       the  fundamental 

1        Jt  differential    for 


differential  formula, 


(21) 


x(t)  =  c1A1  +  c2A2  +  ooeo  +  cnAn  = 


n 


where  Ai(t)  = 


An 

A2i 


A 


ni 


is  the  solution  to  (12)  such  that 


A±j(T)  =  jl,  i  =  j 
,0,  i  J*   J 


^  Xi  =  Jl  °3AiJ  ' 
We  write  (14),  using  (21),  as 

T 


(A-l)      [  X'  6x 


(t)  G  6p  dt 
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From  (20)  and  (21)  we  get  the  m  equations 

(A-2)     A'(t)  G(t)  =  0  . 

It  will  be  recalled  that  (12)  with  (A-2)  constitute  the 
Euler-Lagrange  equations 0 

Additional  relationships  between  the  quantities  dT, 
6Xm,  6p  and  dc.  will  now  be  developed6  See  [4]0 
A  differential  dT  in  T  results  in  the  differential 

dxi(T)  =  6x±(T)  +  x^T  ,  i  =  l,2,..o,n  , 

or 

dxT  =  6xT  +  x  dT 

giving  the  variations 

(A-3)     6x  =  (dx  -  x  dT)T  . 

(A-3)  with  (A-l)  gives 

(A-4)     X«(dx  -  x  dT)T=  U'dx-Xx  dT]T  =/  [x'( t)G( t)6p( t)  ]dt, 

o 
If  x  is  held  constant  in  (A-2),  which  in  expanded  form  is 

i      '  a' 

X'G  =  (0^+  c2il2+  0OO000+  cnAn)  G(x,p,t)  =  0  , 

we  can  write 

(c/jfS  «„  *  A'lG  dCl)  +.....<Cn/\'nf§  6p  +  A'nG  don)  =  0, 
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where  we  define  the  m  x  n  matrix 

§^  =  §^1 •§§>• 

In  terms  of  (21)  this  is  the  set  of  m  equations 

U-5)     X  ~6p  +  A'jG  dc1  +  A'2G  dc2  +  ....+ /A'm  G  dcR  =  0. 

Solving  (A-5)^for  ftp-,,  6p2, ...„„. ,  6p  and  substituting 

v  i 
into  (A-4)  we  get  n  equations,  in  the  n+1  variables 

T,c^, ,c   ,  of  the  form 

/      n  i 

( A-6)     [  Ajdx  =2  Iiidci  +  ^j*  dT  ]T  ,  i  =  1, . . . .  ,n  . 

The  differential  of  the  constraint  h,  is 

(A-7)     dhk  =  J  §§k  dx±  +  §|kdT  L,  ,  k  =  1,2, . . .  „  ,N  . 

For  admissible  curves  h,   is  specified.  The  end  point  must 
therefore  satisfy  the  relation 

(A-8)     dhj^  =  0,  k  =  1,2, ,N. 

In  vector  notation  (A-7)  and  (A-8)  may  be  expressed  by 

(A-9)     [  Vhk  •  dr  ]T  =  0;  k  =  1,2,.. ...,N, 

where  7h,  is  a  gradient  in  x,t  space. 

For  a  maximum  value  of  g(x,t)T  we  similarly  have 

(A- 10)    [  Vg°dr  ]T  =  0 

in  x,t  space  indicating  the  differential  of  the  terminal 
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point  must  be  such  as  to  result  in  displacements  (geomet- 
rically speaking)  normal  to  the  Vg  in  x,t   space,  that  is, 
tangent  to  the  surface  g(x,t)T  =  g ( maximum) . 

Since  the  functions  g(x,t)   and  h(x,t)   are  assumed 
to  be  independent  the  N+l  by  n+l  matrix  of  coefficients 
from  (A-9)  and  (A- 10) 


(A- 11) 


M  = 


S&l      ill!  ».«.  2%  fihjL 


ax,  ax 


fib,  &b 

ax1  ax^ 

&£  l£ 

ax,  ax. 


h 


>x£  at- 


ax£  aF 
5xn  at 


has  rank  N  +  10 

The  transversal  condition,  [4]  and  [2]  page  196,  may 
be  stated  as  the  condition  that  the  rank  of  the  N+2  by 
n+l  matrix 


U- 12) 


MT  = 


M 


Xl  X2 


e  ©  o  o 


n 


X'x 


t  =  T 


be  N+l  and  the  rank  of  the  matrices  obtained  by  omitting 
either  of  the  last  two  rows  also  be  N+l„  Row  number  N+2 
of  M+  is  the  row  of  coefficients  of  dx  and  dT  of  (A-4). 
These  conditions  on  (A-12)  yield  n-N  relations  involving 
the  terminal  values  of  x,  X,  and  t  which  must  be  satis- 
fied by  a  solution  which  gives  at  least  a  stationary  value 
to  g(x,t)To   We  designate  them 
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( A— 13)     H(,x,X,t)m  "■  0 ,    s  —  Ltji  .mi  o.jii-  Wo 

The  differential  relations  derived  from  (A-13),  at  the 
terminal  point ,  are 

U-14) 


dHc  =2  £!s  dx1  +  Z2£5s^i  dC^/Z^s  JU  i5s\    dT 
s        i   dx±        x     i   j   9*i9Ci  i9Xi  6t  ' 

S    —    llOa  g  (   g  o  o  oll'll  j 

where  Xi  is  defined  as  in  (21) 0 

Equations  (A-6),  (A-9)   and  (A-14)  are  2n  equations 
in  the  2n+l  differentials  dx^  dT  and  dC.,  i  and  j  =  l,.,n. 

A, 2     Calculation  procedure 

(a)  A  solution  of  the  problem  is  characterized  by  the  n+1 
parameters  T,C-,,C2, . . .  0C  .     We  arbitrarily  choose  one 
relation  among  the  C.0  We  might  choose  iL(C.)  =  1  or  just 
designate  one,  say  C  =  10  We  must  be  careful  only  that  the 
designated  C.  does  not  turn  out  to  be  negative  or  zeroc  We 
then  guess  a  set  of  values  for  the  remaining  parameters  and 
compute  simultaneously  a  solution  to  the  original  differential 
equations  (3) ,  the  adjoint  system  (12)   and  the  Euler 
equations  (A-2)   or  (20) 0   This  yields  an  extremal  E 
which  does  not,  in  general,  satisfy  (4)  or  (A-13)0 

(b)  Simultaneously  calculate  the  correction  integrals  (A-6). 

(c)  Treat  the  differential  relations  (A-5)   through  (A-13) 


as  differences  evaluated  at  the  end  point  of  E  .  We  then  have 
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(A-15)     Ah.  =  hj  -  h.(t,x)rjn  ,   i  =  1,. ...... N 

AH.  =  -H.(T),   i  =  N+l j . • e • • • 0 • • • jH. 

(d)  Solve  for  the  n  differences  A^, « . .  .,AC  1,  AT  and 
apply  to  the  previous  estimate  for  these  parameters .   If  the 
original  estimate  was  close  enough  the  new  parameters  will 
result  in  a  more  nearly  admissible  extremal,, 

(e)  The  process  is  repeated  until  the  sequence  of  extremals 
converges  to  an  admissible  extremal  within  the  limit  of  some 
convergence  criterion.   Additional  analysis  is  required  to 
distinguish  between  maxima,  minima  and  stationary  values. 
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Appendix  B  THE  FL-METHOD 

Bcl  Outline  of  analysis 

We  have  previously  defined  the  relationships 

(3)  x  =  f(x,p,t)        n  state  variable  equations 

(4)  h(x,t)T  =0         k  constraints 

(10)      6x  =  F  5x  +  G  6p     variational  relationships 

(12)     A=  -F'/\  adjoint  differential  equations 

T 

(14)      [^6x]T=  /  A'g  6p  dt  fundamental  differential 

o         formula 

Following  Bryson[6]  we  define  one  of  the  functions  of 
(4)  as  a  stopping  condition, 

(B-l)     Jl(x,t)T  =  0 

which  produces  a  terminal  value  T  for  t,  thus  eliminat- 
ing the  need  for  guessing  T0 

In  section  (105)  we  developed,  from  solutions  of  the 
adjoint  system, 

(17)      6h.(f)  =  /A'  G  6p  dt,   j  =  1,2, ,k, 

J        o  uj 

or  equivalently 

T 
(B-2)     6h(T)  =  £/(h     G  6p  dt 
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where  Ah(t)   is  defined  as 


Ahct)= 


Qll   n21      nkl 

n12    n22       nk2 


lln 


xkn 


nxk  . 


We  also  developed 


i 
(18)      6g(T)  =  £a*     G  6p  dt  . 

^       Q       & 

We  now  add,  for  (B-l) ,  a  similar  relationship 


(B-3)     «i(T) 


=  X4 


G  5p  dt 


where  AQ  is  the  solution  to  (12)  such  that 


A^(T)  =  VXA(T)  , 


where  v*x  denotes  a  gradient  in  x-space. 

For  small  perterbaticns  the  value  of  T  will  be  changed 
only  a  small  amount  dT,  so  that 


(d-4) 


dg 

=  6g 

+   gdT 

dh 

=  5h 

-i    hdT 

dn 

=  bfi 

+  hdT 

.-/here  g,  h,  A,   are  defined  as 
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g       let      dx1/T 

(B"5)  i  -  ftih  +  ah,\ 

h  ~  (at  +  axxJT 


*-(&*$ 


at      ax/  t 

Accordingly  (18),  (B-2)  and  (B-3)  modified  for  a  change 
dT  in  T  become 


dg  =  AA'g G  6p  dt  +  * 


dT 


(B"6)  /^A'   n 

dh  =  yt  /\.  G  6p  dt  +  h  dT 

o 

T 

dQ  =  ft  An  G  6p  dt  +  n  dT 

o 


In  order  to  limit  the  magnitude  of  the  step  6p 
such  that  the  linearization  (10)  is  reasonable  Bryson  im- 
poses the  limitation 

,T 


(B-7)      (ds)2  =  j£    6p'  6p 


dT 
o 

where  (dS)  is  a  magnitude  arbitrarily  chosen  as  a  reasonable 

upper  limit  to  the  integral  (B-7). 

Throughout  this  paper  a  simplified  version  of  the  FL- 

method,  as  developed  by  Bryson,  has  been  employed  in  that 

variations  at  t  =  t   have  not  been  considered  and  in  that 

o 

his  matrix  of  'weighting  functions',  W(t),  has  been 
neglected. 
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Elimination   of      IT     in    equations   (B-6)    gives 
(B-8)  dg   w  ft  A'ga     C   5p  dt 


c 


r 


(B-9)     dh  =vyC/t'ha  G  6p  dt 


where 


From  (B-7),  (B-8)  ana  (B-9)  we  form  the  linear 

combination 

T 
(B-10)     dg  =  ft   (f{       G  -v'/^  G  -  M6p)  *p  dt  +  v'dh  +,j  (dS)2 


The  second  variation  of  (B-10)  is 

i 
'o 


6Ug)  =  £   (A*   G  -vVt^  G  -  2M6p')  62p  dt 


For  maximum  dg  the  coefficient  of  6  p  =  0  implies 

(B-ll)    6p  =^G-  (/\gn-Ahnv). 

Eliminating     v     and     \x     between   (B-7) ,    (B-9)    and    (B-ll)    we 
have 

(B-12) 


*■                                          /(ds)2-  dh'irJ  dh  - 

6p(t)  =  ±G'(AgcrAh0   Ij£  Ihg)/- f  ♦   G'Ah0   1-1   dh 

^gg"  -"hg  1hh  1hg 


■**   Use  the  plus  sjgn  for  maxima. 
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wnere 


xhh  =  Am  G  a'Am  dt  »    ^g  =  ft%n  G  G*go  dt 


Igg=/tJl'gn°0,^godt 


B.2  Calculation  procedure  for  FL-method 

1.  State  a  reasonable  nominal  control  variable  program, 
p  (t)  which  is  defined  for  the  time  period  (t  ,T)e  In 
many  problems  this  can  be  a  constant  program.  For  ex- 
ample if  the  parameters  are  two  directions,  say  the 
vertical  angle  <()  and  the  bearing  6,  we  could  esti- 
mate p  (t)'=  ^0?G0)»  where  (b  is  any  constant  value 
such  that  0  <  <|)6  <  ir/2  and  6  is  somewhere  in  the  in- 
terval (correct  heading  ±   ir/2) .  Usually  better  estimates 
can  be  made. 

2.  With  p  (t)  solve  system  (3)  determining  T  and  x(t) 
for  tQ  «  t  ,£  T. 

3«,  Solve  the  adjoint  system  (12)  backwards  simulta- 
neously building  Ahn,  Agn,   Ij^,   Ihg,   Igg.   Store 

4nG  ™*  A'mG   • 

4e   On  the  basis  of  x(T)  from  step  two,  evaluate 
h(x,t)T  o   Set  dh  =  -h(x,t)T  . 

5«   Select  a  reasonable  value  for  dS.   If  the  numerator 
of  the  radical  in  (B-12)   is  negative  scale  down  dh 
such  that  the  numerator  is  ^0. 
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6.  Produce  6p*(t)   from  (B-12)   for  t  ^  t  4  1 
simultaneously  forming  p  (t)  =  p  (t)  +  6p  (t)   and 
store  for  the  next  iteration. 

7.  Return  to  step  two  forming  a  new  solution  until  the 
halt  criteria  are  satisfied. 
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Appendix  C   Flow  diagram  for  adjustment  of  (ds) 
Definition  of  symbols 

dh   -the  error  in  meeting  the  constraints  h(x,t)m  =  0. 
e,    -the  magnitude  of  dh  within  which  it  becomes 

reasonable  to  be  concerned  about  extremality. 
e2   -the  convergence  criterion  for  extremality. 
Ai    -a  part  of  (ds)  whose  magnitude  is  adjusted  down  as 

the  solution  approaches  an  extremal. 
D.    -the  denominator  of  the  radical  of  (E-12) , 

dg/ds  in  p- space. 


Calculate  x(T),  h(x,t)T  for  the  nominal 
control  vector  p(t).  Set  dh  =  -h(x,t)  . 


i 


no 

— <■ 


<^dh'dh  <  ^> 


V  yes 


i<  Di-i 


no 


Set  Ai=  0.2  Ai_1 


yes 


V  h-i 


no 


U(dso)2-.(dh'I^  dh+A^)]  <  (^ 


yes 


dh  =  0.5  dh 


y.  no 

set   (ds)2=dh'I^dh+A^ 

)  ^ 

y 

■ 

Calculate  dp   and   return 

"X-   D-  <  62  could  also  be  used  here. 
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